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Physical  and  Chemical  Factors  Affecting 
Contaminant  Hydrology  in  Cold  Environments 

STEVEN  A.  GRANT 


INTRODUCTION 
Definition  of  cold  regions 

Lexically,  the  adjective  “cold”  refers  to  an  environment  in  which  the  ambient  tempera¬ 
ture  is  noticeably  below  body  temperature.  A  cold  region  is  one  in  which  lower  tempera¬ 
tures  have  significant  effects  on  the  natural  environment  or  human  activities.  While  moder¬ 
ate  temperatures  may  be  perceived  as  “cold,”  cold  regions  have  been  typically  defined  by 
the  intervals  in  which  ambient  temperatures  are  below  the  melting  temperature  of  water, 
because  it  is  at  these  temperatures  that  the  effects  of  cold  temperatures  are  most  pronounced. 
Among  the  criteria  that  have  been  used  to  delineate  cold  regions  are  (Bates  and  Bilello 
1966) 

1.  Air  temperatures  below  0°C  or  -18°C  (32°F  and  0°F,  respectively)  that  had  a  50%  like¬ 
lihood  of  being  observed  annually 

2.  Mean  annual  snow  depth 

3.  Ice  cover  on  navigable  rivers 

4.  Isolines  based  on  permanence,  depth,  and  continuity  of  frozen  ground. 

Traditionally,  the  maps  made  from  plotting  isograms  based  on  these  criteria  delineate 
the  changes  in  human  activities  due  to  cold.  For  example,  isograms  based  on  frozen  ground 
delineate  changes  in  the  construction  requirements  for  building  footings.  Isograms  based 
on  ice  cover  speak  to  the  navigability  of  the  waters  during  some  portion  of  the  year. 

This  report  discusses  how  cold  temperatures  affect  contaminant-transport  modeling; 
cold  regions  should  be  delineated  by  the  cold-regions  phenomena  necessary  to  make  the 
model  a  valid  representation  of  the  pertinent  chemical,  physical,  and  microbiological  pro¬ 
cesses  that  determine  the  fate  of  contaminants  in  the  environment.  Cold  temperatures  af¬ 
fect  the  physics  of  contaminant  transport  by  freezing  the  water  in  the  ground  (sometimes  to 
great  depths)  and  by  blanketing  the  ground  (seasonally  or  permanently)  with  snow.  Since 
no  maps  have  been  drawn  delineating  the  cold-regions  effects  on  contaminant  hydrology, 
the  natural  delineations  would  be  those  for  snow  covers  and  extents  of  ground  freezing. 

Extent  of  cold  regions 

In  comparing  the  hydrologic  systems  of  cold  regions  with  those  of  warmer  areas,  two 
factors  distinguish  the  former  systems: 

1.  Some  of  the  annual  precipitation  occurs  as  snow,  which  completes  its  role  in  the  hy¬ 
drologic  cycle  during  a  comparatively  brief  period  as  snowmelt. 

2.  The  ground  freezes  to  some  depth.  Ground  freezing  reduces  the  soil’s  permeability  and 
its  water-storage  capacity.  Accordingly,  ground  freezing  dramatically  decreases  the 
soil’s  infiltration  rate  and  just  as  dramatically  increases  run-off  from  the  soil  of  water 
from  rain  and  snowmelt. 
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Figure  1.  Cold-regions  boundaries  as  determined  by  snow 
depths .  (From  Bates  and  Bilello  1966.) 

To  be  hydrologically  important,  a  cold  region  must  receive  an  appreciable  proportion  of 
its  precipitation  as  snow,  and  its  ground  must  be  frozen  so  that  the  infiltration  of  the  melt¬ 
ing  snow  is  limited.  In  the  northern  hemisphere,  these  two  criteria  can  be  applied  to  delin¬ 
eate  the  extent  of  cold  regions  for  the  purposes  of  contaminant-transport  modeling. 

Snow  cover 

Bates  and  Bilello  (1966)  developed  maps  with  isograms  of  various  annual  snow  depths 
(Fig.  1). 

Much  of  the  United  States  has  average  annual  snow  depths  of  0.3  m  (1 2  in.)  or  more.  This 
includes  all  or  part  of  the  following  states:  Alaska,  California,  Connecticut,  Idaho,  Illinois, 
Indiana,  Massachusetts,  Maine,  Michigan,  Minnesota,  Missouri,  Montana,  Nebraska,  Ne¬ 
vada,  New  Hampshire,  New  Jersey,  New  York,  North  Dakota,  Ohio,  Pennsylvania,  South 
Dakota,  Utah,  Vermont,  Washington,  Wisconsin,  and  Wyoming.  The  isogram  in  Figure  1 
indicates  the  annual  quantity,  but  not  the  annual  fraction,  of  precipitation  that  typically 
arrives  as  snow  and,  therefore,  does  not  indicate  the  relative  importance  of  snow  to  a 
regions  hydrology.  Accordingly,  the  isogram  in  Figure  1  may  not  adequately  delineate 
snow-affected  areas  in  drier  regions. 

Frozen  ground 

Three  major  classes  of  frozen  ground  are  recognized:  seasonally  frozen  ground,  discon¬ 
tinuous  permafrost,  and  permafrost. 

Seasonally  frozen  ground.  The  potential  inadequacy  of  snow-depth  isograms  to  delineate 
cold,  dry  areas  is  borne  out  by  the  map  of  permafrost  and  frost-affected  areas  (Fig.  2),  also 
by  Bates  and  Bilello  (1966).  Understandably,  the  southernmost  isogram,  that  for  frost  pene¬ 
tration,  is  similar  to  that  for  snow  cover.  However,  several  areas  (mostly  in  the  western 
United  States)  with  low  rates  of  annual  precipitation  are  now  included.  These  areas  are  all 
or  parts  of  Arizona,  Colorado,  Kansas,  Maryland,  Missouri,  New  Mexico,  and  Rhode  Is¬ 
land. 

Discontinuous  permafrost  In  the  United  States,  discontinuous  permafrost  is  found  only  in 
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Figure  2.  Cold-regions  boundaries  as  determined  by  frozen 
ground .  (From  Bates  and  Bilello  1966.) 


Alaska.  Discontinuous  permafrost  or  permafrost  is  found  under  most  of  the  land  surface  of 
that  state. 

Permafrost  Areas  of  permafrost,  in  which  no  seasonally  thawed  land  occurs,  are  found  in 
the  northernmost  parts  of  Alaska. 

Importance  to  contaminant-transport  modeling 

A  model  is  a  mathematical  expression  of  the  developer’s  understanding  of  the  system 
being  simulated.  Accordingly,  contaminant-transport  models  valid  for  cold  regions  differ 
from  those  appropriate  for  warmer  climates  because  the  system  being  described  is  conspic¬ 
uously  different.  The  following  modeling  aspects  differentiate  contaminant-transport 
models  valid  in  cold  regions: 

1.  Much  of  what  is  understood  about  the  hydrology  of  cold  regions  is  qualitative.  Be¬ 
cause  there  have  been  comparatively  few  studies  in  cold  regions,  the  physical,  chem¬ 
ical,  and  biological  processes  that  determine  the  contaminant  hydrology  of  cold  re¬ 
gions  are  not  completely  understood.  This  lack  of  quantitative  knowledge  hinders 
the  development  of  physically  based  contaminant-transport  models  valid  for  cold 
regions.  Much  basic  research  is  needed  to  support  the  development  of  simulation 
models  that  are  appropriate  for  these  regions. 

2.  The  energy  flow  and  aqueous-solution  phase  transitions  must  be  included  explicitly 
in  the  development  of  solute  transport  models. 

3.  In  unfrozen  porous  media,  liquid  water  moves  largely  in  response  to  hydrostatic  gra¬ 
dients.  In  frozen  porous  media,  liquid  water  movement  in  response  to  osmotic  or 
thermal  gradients  may  be  dominant. 

4.  The  physical  chemistry  of  electrolyte  solutions  in  natural  porous  media  above  freez¬ 
ing  temperatures  is  well  understood,  and  models  describing  them  are  generally 
available.  For  systems  below  0°C,  there  are  comparatively  few  physical  chemical  data 
(aside  from  freezing-point  depression  determinations)  for  even  simple  electrolyte  so¬ 
lutions. 
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5.  Much  of  the  groundwater  in  areas  with  permafrost  is  brackish  and  impotable.  Un¬ 
lucky  contamination  of  currently  exploited  aquifers  that  are  potable  can  effectively 
deprive  communities  of  potable  subterranean  water  supplies. 

How  these  effects  are  incorporated  in  mathematical  models  of  contaminant  transport 
can  be  outlined  with  the  differential  equations  that  form  the  basis  of  many  of  these  models. 
Water  in  a  porous  medium  flows  in  a  linearly  proportional  response  to  the  hydrostatic 
gradient  across  that  medium.  This  observation  has  been  formalized  as  Darcy’s  law,  the 
vector  form  of  which  is 

v  =  ~  •  (Vp  —  pgVz)  =  -K  •  Vh  (1) 

where  v  =  Darcian  flow  velocity  vector  (m  s_1) 
k  =  permeability  tensor  (m2) 
v  =  viscosity  of  fluid  (Pa  s) 
p  =  fluid  pressure  (Pa) 
p  «  fluid  density  (kg  rrr3) 
g  =  acceleration  due  to  gravity  (m  s-2) 
z  =  vertical  distance  above  datum  (m) 

K  =  hydraulic  conductivity  tensor,  which  in  this  treatment  is  assigned  dimensions 
of  m2  Pa-1  s_1 
h  =  total  head  (Pa). 

(The  formulas  presented  here  are  from  Mangold  and  Tsang  [1991].)  As  will  be  discussed  in 
later  sections,  cold  temperatures  dramatically  affect  three  of  the  parameters  in  eq  1:  v,  p, 
and  K.  The  viscosity  of  the  aqueous  solution  in  freezing  soils  is  affected  directly  by  lower¬ 
ing  temperatures,  which  cause  the  viscosity  of  pure  water  to  increase  by  an  order  of  magni¬ 
tude,  and  by  solutes  excluded  from  ice  forming  in  the  soil  solution,  which  cause  the  viscos¬ 
ity  of  the  remaining  liquid-water  solution  to  be  increased  still  further.  The  pressure  of  vici¬ 
nal  water  in  frozen  porous  media  is  controversial  and  an  active  area  of  basic  research 
(Zheng  et  al.  1991).  It  is  clear  that  the  liquid-pressure  gradients  accompanying  differences 
in  temperature  may  be  much  more  important  than  pressure  gradients  due  to  elevation  or 
water  content  (see,  for  example.  Perfect  et  al.  [1991]).  Finally  as  ice  occupies  some  of  the 
pore  space,  the  conductivity  of  the  porous  material  is  reduced. 

Once  the  flow  of  water  is  established,  the  more  challenging  problem  of  solute  flows 
through  porous  media  can  be  addressed.  In  many  cases,  the  appropriate  partial  differential 
equation  to  describe  the  transport  of  the  solute  j  is 

v#(vCj)-v,(d#Vcj)  =  cJ  +  Qci  (2) 

where  Cj  -  aqueous  concentration  of  jth  solute  (mol  rrr3) 

D  =  dispersivity  tensor  (m2  s'1) 

0  =  porosity  of  porous  medium  (dimensionless) 

Qcj  =  source-sink  term  for  jth  solute  (mol  nr3  s_1). 

As  would  be  expected,  most  of  the  geochemical  reactions  that  contribute  to  the  terms  Qcj 
are  affected  by  temperature.  The  most  pronounced  of  these  effects  are  on  the  chemical- 
thermodynamic  state  of  the  solutes  and  consequently  on  their  solubility,  miscibility,  speci- 
ation,  and  reactivity. 

AQUATIC  CHEMISTRY  BELOW  0°C 

The  transport  of  solutes  and  non-aqueous-phase  liquids  in  any  porous  medium,  frozen 
or  unfrozen,  is  directly  affected  by  the  physical-chemical  properties  of  the  solvents  and 
solutes  in  the  system.  These  properties  can  be  measured,  but  in  many  cases  the  careful 
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application  of  chemical-thermodynamic  information  can  spare  the  modeler  from  making 
many  tedious  measurements  in  the  laboratory.  While  few  systems  have  been  studied  as 
extensively  as  aqueous  systems,  there  is  a  paucity  of  data  at  subzero  temperatures.  This 
section  discusses  the  physical-chemical  properties  of  the  solvent,  water,  as  a  pure  liquid. 
After  a  brief  discussion  of  the  relevant  chemical  thermodynamics,  the  section  concludes 
with  a  discussion  of  the  properties  of  aqueous  solutions  at  subzero  temperatures. 

Physical  properties  of  liquid  water 

The  state  and  movement  of  water  in  frozen  and  unfrozen  ground  is  affected  by  several  of 
its  physical  properties,  all  of  which  change  with  temperature  to  some  degree.  Some  of  these 
temperature  effects  are  trivial,  but  others  must  be  taken  into  account  explicitly  when  mod¬ 
eling  solute  transport. 

The  most  up-to-date  tabulations  of  the  physical-chemical  properties  of  liquid  water 
have  been  presented  in  the  Journal  of  Physical  and  Chemical  Reference  Data  (a  recent  cu¬ 
mulative  index  is  presented  on  page  1368  of  volume  20).  Many  of  these  properties  can  be 
modeled  conveniently  with  the  numerical  equation-of-state  package  developed  by  Hill 
(1990). 

Most  of  the  physical  properties  of  supercooled  water  have  been  measured  to  roughly 
-35°C.  The  temperature  trends  of  these  physical  properties  are  consistent  with  a  limiting 
behavior  for  many  physical  properties  of  liquid  water  at  a  unique  temperature,  apparently 
-46°C.  Speedy  (1987)  fitted  three  properties  of  supercooled  water— cubic  expansion  coeffi¬ 
cient  (a),  isothermal  compressibility  (1C7),  and  heat  capacity  at  constant  pressure  (Cp) — to  an 
empirical  relation: 

X=  £  B e"  +2 Cxt~m ,  (3) 

n=0 

where  X  is  the  physical  property  being  described.  The  empirical  coefficients  Ex  and  Cx 
have  been  fitted  to  the  model  from  experimental  data.  The  parameter  e  is  a  reduced  temper¬ 
ature  [=  ( T  -  Ts)/Ts],  Ts  being  a  limiting  temperature  that  is  assumed  to  be  exactly  227.15  K 
(-46°C).  The  fitted  parameters  of  the  Speedy  model  are  presented  in  Table  1. 


Table  1.  Coefficients  to  Speedy’s  (1987)  empirical  equation  for  cal¬ 
culating  heat  capacity,  thermal  expansivity,  compressibility,  and 
viscosity  of  supercooled  water. 


Parameter 

Cubic  expansion 
coefficient 
ax  103 
(K) 

Isothermal 

compressibility 

kt 

(bar1) 

Heat  capacity  at 
constant  pressure 

cp 

OK-1  moH) 

Viscosity 

V 

(Pa  s-1) 

Cx 

-0.80 

20.0 

14.2 

26.312 

40) 

1.802  180  3 

4.120 

24.952 

-144.565 

Bx 

-0.941  698  0 

-1.130 

128.281 

1,239.075 

Bx 

0.905  507  0 

77.817 

14.2 

-8,352.579 

43) 

-0.80 

-78.143 

-221.405 

31,430.760 

44) 

54.290 

-64.812 

48,576.798 

Max  resid 

1.2  ppm 

0.2% 

0.03% 

1.9% 

Density 

Water  at  equilibrium 

Compared  with  its  other  physical-chemical  properties,  the  density  of  liquid  water 
changes  little.  Under  a  pressure  of  0.1  MPa,  the  density  of  air-free  water  is  at  its  maximum 
at  3.98°C.  Under  these  conditions,  the  density  of  water  is  1000  kg  m~3.  This  is  the  density  of 
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aqueous  solutions,  which  is  most  often  assumed  in  contaminant-transport  models,  even 
though  water  at  25°C  under  atmospheric  pressure  has  a  density  of  997.02  kg  nr3,  a  relative 
change  of  2.98%  from  its  maximum.  At  the  triple  point  of  water  (0.01°C),  liquid  water  has  a 
density  of  999.78  kg  m-3,  a  relative  change  of  -0.022%  from  its  maximum. 


Supercooled  water 

Speedy  (1987)  has  presented  equations  that  have  been  fitted  to  the  thermodynamic 
properties  of  supercooled  water.  His  equation  for  the  density  of  water  is 


p  =  psexp\-Ts 


f5 


I^B«n)en+,+2Cael/2 


n=0 


(4) 


where  ps  is  a  reference  density  and  both  and  Ca  are  empirical  parameters,  the  values  of 

which  are  presented  in  Table  1 .  A  global  fit  of  eq  4  to  the  density  data  of  Hare  and  Sorensen 


(1987)  yielded 

the  following  parameter  estimates: 

Ps 

0.901  532  859  3 

Qt 

-0.001  176  165  2 

40) 

0.003  844  238  2 

4° 

-0.015  727  076  1 

42) 

0.074  406  461  4 

43) 

-0.140  643  265  3 

The  largest  deviation  of  the  fitted  curve  from  measured  data  is  0.20%. 


Viscosity 

The  viscosity  of  liquid  water  increases  exponentially  with  decreasing  temperature. 
While  the  International  Association  for  the  Properties  of  Water  and  Steam  (IAPWS)  has 
released  an  elaborate  expression  for  calculating  the  viscosity  of  water  as  a  function  of  pres¬ 
sure  and  temperature,  a  much  simpler  expression  may  be  used  accurately  for  common 
pressures  and  temperatures  between  -40  and  25°C.  Experimental  viscosity  data  have  been 
fitted  to  the  empirical  Vogel-Tammann-Fulcher  (VTF)  equation: 


v  =  v0  exp 


r  B 
J  +  t0 


(5) 


where  v0  =  0.028  556  Pa  s 
B  =  509.53°C 
t  -  temperature  (°C) 

t$  =  empirical  constant  equal  to  123.15°C. 

A  modification  of  this  equation  is  able  to  fit  the  available  data  more  closely: 


v  =  v0  exp 


(- 

v+l 


B-> 


+*o  (?+*or 


h. 

(*+/o)V 


(6) 


where  vq  =  0.000  046  01  Pa  s 
Bx  =  3068.6°C 
Bz  =  -3.3775  xl05oC2 
JS3  =  1.4781  x  107  °C3. 

This  relation  is  plotted  in  Figure  3,  showing  its  agreement  with  measured  values. 

In  the  temperature  range  of  -35°C  to  0°C,  the  viscosity  of  supercooled  water  may  be 
fitted  to  Speedy’s  (1987)  limiting-temperature  empirical  equation  (eq  3).  This  relation  is 
presented  in  Figure  4. 
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Figure  3.  Supercooled  water  viscosity  modeled  by  modified 
Vogel-Tammann-Fulcher  (VTF)  equation . 


Temperature  (°C) 

Figure  4.  Supercooled  water  viscosity  estimated  by  Speedy's  (1987) 
empirical  formula. 


In  addition  to  temperature,  the  viscosity  of  aqueous  solutions  is  affected  by  the  type  and 
concentration  of  solutes.  These  effects  are  discussed  in  standard  references  (e.g.,  Harned 
and  Owen  1958,  Robinson  and  Stokes  1959).  Vand  (1948)  derived  an  equation 


lavrei  = 


A^C 

1  -Qe 


(7) 


where  vrej  is  the  ratio  of  the  solution’s  viscosity  to  that  of  the  water. 

This  author  found  that  an  extension  of  the  Vand  equation  gave  improved  fits  to  viscosi¬ 
ties  of  NaCl  solutions: 


{Aq+AxT)c 

nVrel  i-{Qb+Q{Ty 


(8) 


This  equation  was  fitted  to  the  data  of  Lyle  and  Hosking  (1902),  who  reported  the  viscosity  of 
NaCl  solutions  with  initial  concentrations  (at  15°C)  of  0.1,  0.2, 0.5, 1.0,  2.0,  and  4.0  mol  dm-3. 
The  following  parameter  estimates  (and  standard  errors  of  the  estimates)  were  obtained: 
Ao  =  -0.116  90  (0.003  97),  Aj  =  0.000  579  (0.000  012),  Qo  =  0.440  930  (0.003  467),  and  = 
-0.001  065  909  6  (0.000  000  000  00). 
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Figure  5.  Predicted  viscosity  of  aqueous  NaCl  solutions  as  affected  by  temperature  and 
solution  concentration. 


This  relation  and  an  equation  that  describes  the  viscosity  of  supercooled  water  (e.g.,  eq  6 
or  eq  3)  allow  one  to  predict  the  viscosity  of  electrolyte  solutions  at  low  temperatures, 
though  predictions  at  high  molalites  or  low  temperatures  may  be  subject  to  significant 
extrapolation  errors.  To  the  author’s  knowledge,  there  are  no  measurements  of  viscosities 
of  electrolyte  solutions  at  subzero  temperatures  in  the  published  literature.  A  plot  of  the 
likely  interaction  between  temperature  and  electrolyte  molality  on  the  viscosity  of  aqueous 
solutions  is  presented  in  Figure  5. 

Interfacial  tension  of  liquid  water  against  its  vapor 

The  interfacial  tension  of  water  against  air  (in  N  nr1)  from  the  triple  point  of  water 
(0.01°C)  to  the  critical  point  of  water  (374.15°C )  has  been  fitted  to 

in  which  the  parameters  have  the  following  recommended  values:  yo  =  0.2358  N  nr1,  T\  = 
647.15  K,  u  =  1.256,  and  v=  0.625.  While  there  have  been  few  reported  measurements  of  ywa 
below  0°C,  these  data  indicate  that  eq  9  is  approximately  valid  to  at  least  -8°C  (Haar  et  al. 
1984).  The  complexity  of  eq  9  belies  the  fact  that  the  interfacial  tension  of  water  against  its 
vapor  is  virtually  linear  from  -10°C  to  50°C,  as  can  be  seen  in  Figure  6, 

Interfacial  tension  of  liquid  water  against  ice 

Apparently,  the  most  precise  measurement  of  the  interfacial  tension  of  liquid  water 
against  ice  is  that  of  Ketcham  and  Hobbs  (1966),  who  reported  a  value  of  y fW  (p  =  0.101325 
MPa,  t  =  0°C)  =  0.033  N  nr1.  A  compilation  of  other,  earlier  measurements  was  presented 
by  Jellinek  (1972). 

Heat  capacity 

Heat  capacity  is  the  fundamental  physical-chemical  quantity  by  which  changes  in  en¬ 
tropy  and  enthalpy  with  temperature  can  be  calculated.  Given  the  chemical  and  biological 
importance  of  water,  it  is  understandable  that  equilibrium  heat  capacities  of  pure  liquid 
water  and  ice  phases  have  been  measured  precisely  for  most  temperatures  and  pressures. 
In  recent  years,  the  heat  capacities  of  supercooled  water  have  been  measured  with  ever- 
increasing  precision. 
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Figure  6.  Measured  values  of  interfacial  tension  of  water  against  its 
vapor  at  a  range  of  temperatures. 


Water  at  equilibrium 

The  heat  capacity  of  liquid  water  can  be  calculated  accurately  from  the  equation  of  state 
for  liquid  water.  The  equation-of-state  models  for  water  are  comprehensive  and  algebra¬ 
ically  complex.  Fortunately  there  are  several  computer  programs  available  (Harr  et  al. 
1976,  Hill  1990)  that  can  be  used  to  calculate  accurately  the  thermodynamic  properties, 
including  heat  capacity. 

Supercooled  water 

As  noted  above,  Speedy  (1987)  has  presented  equations  for  the  various  thermodynamic 
properties  of  supercooled  water.  His  expression  for  the  constant-pressure  heat  capacity  of 
supercooled  water  is 

Cp=  I  Bg*e®  +  20:  6-1/2.  (10) 

n=0  p  P 

The  regression  estimates  of  these  parameters  are  presented  in  Table  1.  Figure  7  presents  the 
isobaric  heat  capacity  of  water  as  calculated  by  the  models  of  Hill  (1990)  and  Speedy  (1987). 


Figure  7.  Heat  capacity  of  water  as  predicted  by  Hills  (1990)  equa¬ 
tion-of-state  model  (for  t<0°C)  and  Speedy's  (1987)  empirical  for¬ 
mula  (for  t  <  0°C). 
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Chemical  properties 
Thermodynamic  properties 

To  calculate  the  chemical  thermodynamic  properties  of  aqueous  solutions  below  0°C,  the 
relevant  properties  of  the  solvent  must  be  calculated.  The  conventional  reference  tempera¬ 
ture  and  pressure  (298.15  K  and  0.1  MPa)  are  inappropriate  because  of  Speedy’s  (1987)  em¬ 
pirical  model  upper  limit  of  applicability  of  273. 1 5  K  and  were  determined  under  a  pressure 
of  1  atm.  Initially  then,  the  chemical  potential  of  water  at  273.15  K  and  0.101325  MPa  must 
be  calculated. 

Chemical  potential  of  water  at  273.15  K.  To  state  the  chemical  potential  of  water  at  0°C,  one 
needs  to  calculate: 


ta20(l)(^TW’Prw)  M'H20(l)(^rc’Prc) 

where  Trc  =  298.15  K 
T™  =  273.15  K 
Ts  =  227.15  K 
Prc  —  0.1 

p™  =  0.101325  MPa. 

Equation  11  can  be  evaluated  with 

(^rw  *  Prw )  “  HB (^rc » Prc ) 


T 

1rw 

=  -^(Trc,prc)(Trw-Trc)+  J  CpRdT  (12) 

Tc 

^rw  P*  Prw  . 

-Sw  J  -^dT+  J  Vfe  (Trw,p)dp 

^rc  Prc 

where  j^b  Is  the  chemical  potential  of  pure  B  Q  mol-1),  Sb  is  the  molar  entropy  of  pure  B 
Q  K-1  mol-1),  and  v%  is  the  molar  volume  of  pure  B  (m3  mol-1). 


Equation  12  is  equivalent  to 

M-B^rw  *  Prw)  ”  Hb(^tc  *  Prc)  =  (^rc » PrcX^rw  ”  ^rc)  +  j^B^rw  » Prc)  “  ^b(^tc  * Prc)J 


SB*(Trw>prc)-^(Trc,prc)]+  j  v;(Trv/lP)dp 

Prc 


(13) 


where  Hb  is  the  molar  enthalpy  of  pure  B  fl  mol  l). 

Standard  reference  tables  (e.g.,  Chase  et  al.  1985)  provide  some  of  the  required  values: 

^H20(l)(^rc>Prc) =  69.950  JK_1  mol-1. 

The  equation-of-state  model  for  water  (Hill  1990)  can  be  used  to  calculate  some  of  the  terms 
directly: 

SH2O(l)(rrw.Prc)-^I2O(1)(rrc.Prc)  =  -6-616  04  J  K1  mol-1 
WH20(l)(^rw-Prc)_  WH20(l)(^rc>Prc)  =  _1-88547  kJ  mo1  ' 

^H20(l)(^rw.Prc)  =  0.018  01809  dm3  mol-1 
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The  most  recent  accurate  measurement  of  the  enthalpy  of  fusion  for  water  to  ice  appears  to 
be  that  of  Haida  et  al.  (1974),  who  reported 


A'sH^ofT™, ftc)  =  6006.8  J  mor1. 
By  the  relation  (McGlashan  1979) 

AaPHB*  =  ra+pA|<£ 


we  find  that 

21.99085  J  K_I  moP1 
which  allows  the  evaluation  of  eq  12: 


M'B(^rw>/?rw)  M'B^C’A'c)  —  (^ic»  ArcX^rw  ^rc)^B  (^rw’Prc)  ^b(^tc’ A-c) 

Prw  Prw  _J 

-7^Wre)-£frc.Ac)]  +  J  vS{Tm,p)dp+  J  FB(7’rw>/,)dp  =  1.66737kJmol 

Ptc  Ptc 

Supercooled  water 

The  standard-state  chemical  potential  of  supercooled  water  at  temperatures  between 
0°C  and  -46°C  under  pressures  of  1  atm  can  also  be  calculated  by  eq  12.  In  this  temperature 
range,  the  changes  in  Cp  ,h20(1)  are  described  by  eq  10.  The  definite  integrals  needed  to 
make  these  calculations  are  as  follows: 


?C!„, oa,..  W-tt)  [Bg-3Bg»,6Bg]lM^ 

J  -^r~~dT-  47?  +  37? 


27? 


»-2Bg+3Bg)-4BW](Tf-rrw) 


+  2Cc  arctan 


) 


-2  arctan 


--1 


B®  -BjP  +  B<?>  -B<?>  +B^]x[ln(Tf)-ln(Trw)]  (17) 


and 


Tf 

I 

Tr 


f  r*  r\T —( d(1)  .  d(2)  d(3)  ,  p(4) 
J  Cp.H20(l)  d  1  -  ^C„  "  BC„  +  BC„  -  BCn  +  "c, 


v  . 


51? 


[b»;  -&)  [b®  -3B®  +6B»;](lf'-^) 

+  4Z?  +  3!? 


B$  - 

+L£^_ 


2  B®  +3B®  -4B! 
cp  lp 


.(4) 


!_  j2  \ 

Jrw  J 


2T 


+  2Cr 


(18) 
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Chemical  thermodynamics  of  mixtures  and  aqueous  solutions 

Before  considering  the  properties  of  electrolyte  solutions  at  subzero  temperatures,  this 
discussion  surveys  briefly  the  chemical  thermodynamics  of  mixtures  and  solutions.  Con¬ 
sider  a  mixture  of  N  components: 


nj  ~  n A  +  WB  •••**" nN 

where  np,,  ng . . .  are  the  amounts  of  components  A,  B  . . .  (mol)  and  nj  total  number  of  mole¬ 
cules  in  the  mixture  (mol).  If  the  number  of  molecules  of  one  component  is  much  greater 
than  all  the  other  components,  the  mixture  is  called  a  solution.  The  dominant  mixture  com¬ 
ponent  is  called  the  solvent,  and  the  other  mixture  components  are  called  solutes.  If  the 
mixture  being  discussed  is  a  solution,  by  convention  the  letter  A  is  reserved  to  represent  the 
solvent. 

Consider  now  a  binary  mixture  or  solution: 


nT  =  n  A+WB*  (20) 

For  all  extensive  thermodynamic  properties  of  the  mixture  or  solution,  the  partial  molar 
quantities  can  be  defined;  for  example,  the  partial  molar  volume  of  B  (m3  mol-1)  is 


dVA 

dnB  Jp,T,ni3i* B 


(21) 


where  Vy  is  the  total  volume  of  the  mixture  (m3)  and  3  is  a  logical  symbol  for  “such  that.” 
The  total  volume  of  the  binary  mixture  is  therefore 


Vj  =  nA^A  +  • 


(22) 


We  can  perform  a  similar  operation  on  the  total  Gibbs  energy1  (G  j,  J)  of  the  mixture, 


Gy  =  npfik  + 


(24) 


where  the  partial  molar  Gibbs  energy, 

(25) 

V™B  / p, T.njSfeB 

is  defined  similarly  to  before.  And,  as  before,  the  total  Gibbs  energy  can  be  represented  in 
terms  of  partial  molar  quantities: 


Gy  =  +  %Gb  • 

The  partial  molar  Gibbs  energy  is  also  known  as  the  chemical  potential 

HB=fGB.  (27) 

In  a  mixture,  the  chemical  potential  of  each  component  differs  from  some  defined  standard 
state.  This  difference  defines  the  activity  of  the  mixture  component, 


1  This  quantity  is  defined  by: 

def 

Gj  ~  Hj  —  TSj 

where  HT  is  total  enthalpy  (J)  and  Sj  is  total  entropy  (J  K“!) . 


(23) 
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(28) 


Hb=Hb  +RrinaB, 


where  |if  is  the  standard-state  chemical  potential  of  component  B  and  is  its  activity.  In 
many  cases,  the  activity  of  a  mixture  component  can  be  approximated  by  its  mole  fraction, 
defined  by 


(29) 


or  its  molality,  defined  by 


™  ___5 B 
B  («aMa) 


(30) 


where  Ma  is  the  molar  mass  of  solvent  A  (kg  moH). 

These  approximations  are  imperfect,  so  the  deviation  of  the  mole  fraction  or  molality 
from  the  activity  defines  an  activity  coefficient  (fB  or  Yb),  termed  rational  if  the  reference 
function  is  the  mole  fraction: 


pB=nf  +  J?rin(fBXB),  (31) 

or  molality-based  if  the  reference  function  is  solution  molality: 

pB=pf +  J?rin^»-j  (32) 

where  mjf  is  the  standard  state  molality,  1  mol  kg-1. 

Electrolyte  solutions 

Geochemical  solutions  are  electrolyte  solutions.  Because  electrolytes  dissolve  into 
charged  ions  in  solution,  the  chemical-thermodynamic  description  of  these  solutions  is 
much  more  difficult  than  for  solutions  and  mixtures  of  electrically  neutral  components. 

Consider  the  simple  case  of  an  unhydrated  crystalline  J3VbCv  (cr)  composed  of  vc 
moles  of  anion  CZq  and  BZb  (zb  and  zq  are  the  charge  numbers  of  the  subscripted  ions) 
that  dissolves  completely  in  water  to  form  an  anionic  and  a  cationic  species: 


Bv  Cy  (cr)  =  vB  BZb  (aq)  +  vc  Cz°  (aq). 

B  C 

(33) 

The  single-ion  activity  coefficient  Yb  can  be  defined  by 

-  defYBmB 

*<W  =  me 

(34) 

and  related  to  the  chemical  potential  of  the  single-ion  species  by 

.B,aq)=,t(aq)+^ln(^B). 

(35) 

For  niBC  m°les  of  salt,  which  dissociate  into  v%m% c  moles  of  cation  and  vc/hbc  mo^es 
anion,  the  chemical  potentials  for  the  salt  jibc  in  solution  is  related  to  the  chemical  poten¬ 
tials  for  the  single-ion  species  by 

mBC  f-’-BC(aq)  =  mB  M-B(aq)  +  mC  ^C(aq) 

(36) 

(From  eq  37,  the  salt  6VbCVc  will  be  represented  simply  as  BC).  The  chemical  potentials  of 
the  electrolyte  and  the  single-ion  species  are  related  by 

M'BC  =vbIib+vc1xc 

(37) 
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which  leads  to 


Hbc  =  vB|if  +vC|x§  +vBl?Tln(/nBYB)  +  vci?Tln(mcYc) 

(38) 

and 

M-bc  =  nfc  +(vb  +vc)i?TlnJjnBC(v^B  v^c)1/(VB+Vc)y±BCJ. 

(39) 

By 

^BC  =  vCliC  +  vB^B 

(40) 

the  mean-ionic  activity  coefficient  is  therefore  defined  as 

1 

Y±bc  =f(Ys  Yc  )(vb+vc>  .  (41) 


Figure  8.  Standard-state  chemical  potential  of  NaCl(aq)  as  function  of 
temperature. 


Figure  9.  Mean-ionic  activity  coefficient  of  NaCl(aq)  (p  =  0.1  MPa, 
m  =  1.0  mol  kg-1)  as  function  of  temperature. 
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Chemical  equilibrium  modeling  considerations 

To  model  the  chemistry  of  any  aqueous  electrolyte  solution,  the  chemical  potential  of 
each  modeled  constituent  must  be  estimated  either  explicitly  or  implicitly  The  most  effi¬ 
cient  means  of  estimating  the  single-ion  chemical  potential  of  a  solute  is  by 

^B(aq)  =  Ml (aq)  +  •  (35) 

Both  the  standard-state  chemical  potential  of  the  solute  [  ^f(aq)]  and  the  single-ion  activity 
coefficient  ye  are  affected  by  temperature.  These  effects  are  presented  in  Figures  8  and  9. 
Further,  ice  is  nearly  a  pure  phase.  As  aqueous  solutions  are  cooled  below  their  freezing 
points,  ice  forms,  excluding  solutes  to  the  remaining  liquid-water  solutions.  Accordingly, 
the  aqueous  solutions  at  subzero  temperatures  are  typically  brines,  the  activity  coefficients 
for  which  must  be  estimated  by  elaborate  models  such  as  that  of  Pitzer  and  his  coworkers. 
(For  a  further  discussion,  see  Pitzer  1991  and  Spencer  et  al.  1990.)  Figure  10  shows  the  iso¬ 
thermal,  isobaric  changes  in  the  activity  coefficient  of  a  salt  (NaCl)  with  its  concentration. 
Colligative  properties :  freezing  point  depression 


NaCl  Molality  {mol  kg-1) 

Figure  10.  Mean-ionic  activity  coefficient  of  NaCl(aq)  (p  -  0.1 
MFa,  T  =  298.15  K)  as  function  of  molality. 

While  liquid  water  is  an  almost  universal  solvent,  ice  is  generally  a  very  pure  phase.  The 
freezing  point  of  an  aqueous  solution  is  depressed  proportionately  to  the  concentrations  of 
its  solutes.  The  thermodynamics  of  this  behavior  is  fairly  straightforward.  The  discussion 
here  follows  closely  that  of  Atkins  (1990). 

At  equilibrium,  the  chemical  potentials  of  the  solvent  in  the  liquid  and  solids  phases  are 
equal: 


Ha(s)  =  M-a(1)  • 


(36) 


Since  the  standard  state  for  the  solvent  is  pure  solvent  substance: 


^(DUVAd) 


(37) 


by  eq  29  the  following  relation  can  be  derived: 

M-a  (s)  =  M-a  (1)  +  i?7's+1  In  aA  (38) 

is  obtained,  7B+1  representing  the  freezing  temperature.  If  it  is  assumed  that  aA  =  xA  (i.e., 

4-D. 
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Pa  (s)  =  Pa  (1)  +  RTs+l  lnxA. 


(39) 


As  with  total  thermodynamic  properties  (eq  25),  the  Gibbs  energies,  enthalpies,  and  entro- 
pies  of  fusion  are  related  by 


AlsGH2o(T,p)  =  a\Hh2o(T,p)  -  TA‘sSh2o(T,p).  (40) 

From  this  expression,  we  can  obtain  directly 

h  Pa (s)  -  Pa 0)  _  ^M2Q (T  p)  A's^2o(r,  p)  (41) 

A  RTs+l  RTs+\  R 

Recalling  that  for  a  binary  solution 
*A  +  *B  *  1 

one  can  calculate  the  molar  Gibbs  energy  of  fusion  for  two  systems:  a  pure  solvent  and  a 
solution. 

For  the  pure  solvent  system, 


Ini  Q  Pa(s)~PaO)  ashh2o  V*.P)  4^h2o  (T*.p) 
gjs+l*  RTs+l *  R 

and  for  the  solution  system, 

Infi  --  )  Pa(s)~PaQ)  A‘s^H2q(T'P)  As^H2p(r,P)  (43) 

ln(i  -  xg)  -  Rrp+l  -  RTS+l  R 

where  Ts+{*  is  the  freezing  point  of  the  pure  solvent  system.  Subtracting  one  from  the 
other,  we  obtain 


a‘sHh2o(T,p)  a'sHh2o  (T* ,  p) 
ln(l  xB)  -  RTs+i* 


(44) 


With  small  error,  it  can  be  assumed  that  the  enthalpy  of  fusion  changes  insignificantly  with 
small  changes  in  temperature,  that  is. 


As^HjO  (T,  p)  ~  A*sHh2o  (T* ,  p). 
This  leads  immediately  to 


(45) 


ln(l  -  xB) 


A[sH^0(T,p)f  1 _ 1) 

R  [js+l  js+l*) 


(46) 


If  it  is  further  assumed  that  ln(l  -  xg)  =  -xg  and  that  Ts+1 «  P+t* . 

An  approximate  expression  for  the  freezing-point  depression  is  therefore  provided: 


AT = 


Rrz 


AXoff.p) 


*B 


(47) 


where  AT  =  T8-1"1’*  —  T8+1. 
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Equation  47  suggests  that  freezing-point  depression  is  a  solvent  property.  This  observa¬ 
tion  has  been  formalized  in  the  definition  of  a  cryoscopic  constant,  Kf,  for  a  given  solvent: 

AT  =  Kf /nB  (48) 

which  for  water  is  1.86  K  kg  mol-1. 

The  freezing  point  is  actually  a  phase-equilibrium  line  along  which,  as  the  temperature  of 
a  system  decreases,  a  largely  pure  ice  phase  is  in  equilibrium  with  increasingly  concen¬ 
trated  aqueous  solutions.  Examples  of  this  behavior  can  be  seen  in  Figures  11  and  12,  which 


Figure  11.  Phase  diagram  for  NaCl-fyO  system .  (Courtesy  of  Dr.  G.M.  Marion,  USACRREL.) 


Figure  12.  Phase  diagram  for  CaCl 
H20  system.  (Courtesy  of  Dr.  G.M. 
Marion,  USACRREL.) 
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present  phase  diagrams  of  the  NaCl-H20  and  CaCl2-H20  systems  under  atmospheric 
pressure.  It  is  helpful  to  recall  the  phase  rule  when  examining  these  phase  diagrams.  For 
nonreacting  systems,  the  phase  rule  states 

F-C+2-P  (55) 


where  F  is  the  degrees  of  freedom,  C  is  the  number  of  components,  and  Pis  the  number  of 
phases  (McGlashan  1979) .  A  two-phase  system  (i.e.,  ice-NaCl  solution)  that  is  composed  of 
two  components  (i.e.,  NaCl  and  H2O)  has  (2  +  2-2)  two  degrees  of  freedom,  that  is,  fixing 
two  independent  variables  will  specify  the  state  of  the  system.  Three  independent  variables 
are  most  natural:  pressure,  temperature,  and  composition.  Since  pressure  is  fixed  by  the 
atmosphere,  the  choice  of  temperature  or  composition  will  specify  the  other  independent 
variable.  This  dependence  is  represented  by  the  ice-solution  coexistence  lines  in  Figures  11 
and  12.  Simple  NaCl  solutions  with  molalities  above  5.2  that  are  cooled  below  0°C  will  not 
coexist  with  ice  but  rather  with  a  solid  phase  composed  of  hydrohalite  [NaCl  •  2H2(cr)]. 
The  ice-solution  hydrohalite-solution  coexistence  curves  converge  at  the  eutectic  point, 


the  unique  temperature- 
composition  combination 
where  these  three  phases 
coexist,  a  triple  point  that  is 
specified  entirely  by  the 
number  of  phases  and  the 
lone  independent  variable, 
pressure. 

The  phase  rule  can  be 
extended  to  the  interpreta¬ 
tions  of  the  freezing  behav¬ 
ior  of  natural  waters. 
When  natural  waters  are 
frozen,  solute  effects  must 
be  considered.  In  particu- 


0  200  400  600  800  1000 


Water  (g) 


lar,  several  salts  precipitate  Figure  13.  Phase  diagram  for  seawater.  (Courtesy  of  Dr.  G.M. 
from  freezing  solutions  as  Marion,  USACRREL.) 
the  remaining  liquid  water 

solutions  concentrate.  This  progression  has  been  studied  most  closely  in  the  freezing  be¬ 
havior  of  seawater,  the  properties  of  which  may  be  used  as  an  approximate  model  for  soil 
solutions  in  freezing  and  frozen  ground.  A  phase  diagram  for  seawater  is  presented  in  Fig¬ 


ure  13. 


PHYSICAL-CHEMICAL  PROPERTIES  OF  ICE 
Heat  capacity 

As  for  many  solids,  the  heat  capacity  of  ice  as  a  function  of  temperature  can  be  fitted  to 
the  equation  of  Maier  and  Kelly  (1932): 

Cp,H20(cr.D  =  a  +  br+^-  (50) 

where  a  -  -10.6644  (1.5999)  J  K’1  mol'1 
=  0.1698  (0.0046)  JK-2  moH 
c=  198  148  (28  230)  JK  mol'1. 

The  changes  in  molar  enthalpy  and  entropy  can  thereby  be  calculated  via  the  following 
integrals: 
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Figure  14.  Molar  volume  of  hexagonal  ice  under  1  atm  pressure  at  temper¬ 
atures  near  its  melting  point. 


and 


J  Cp'Hy(cr,1)  dT=a[ln(Tf)-ln(Trw)]  +  6(Tf-Trw)-| 

Trw 


(52) 


Molar  volume 

Water  is  one  of  the  few  substances  for  which  one  of  its  solid  phases  is  less  dense  than  its 
liquid.  The  molar  volumes  of  hexagonal  ice,  the  phase  most  commonly  encountered  in 
nature,  at  temperatures  near  its  melting  point  are  presented  in  Figure  14. 


PHYSICS  OF  FROZEN  POROUS  MEDIA 

Moderately  cool  temperatures  affect  the  physical-chemical  state  of  water  and  its  pro¬ 
pensity  to  flow  through  porous  media.  For  contaminant  transport,  changing  temperatures 
affect  die  chemical  potentials  of  the  solvents  and  solutes  in  the  pore  solutions  and  therefore  the 
solubility  and  reactivity  of  the  solutes.  These  changes  in  the  physical  properties  of  the  sol¬ 
vent  and  chemical  properties  of  the  solutes  can  affect  dramatically  the  mobility  of  contam¬ 
inants  in  cold  regions. 

At  cold  temperatures,  at  which  surficial,  vadose,  and  aquifer  waters  freeze,  the  nature  of 
flow  changes  dramatically.  First,  at  a  macroscopic  level,  because  the  effective  volume  of 
voids  is  reduced,  ground  has  a  much  lower  permeability  when  frozen.  This  affects  the  ba¬ 
sin-scale  hydrology  in  cold  regions  in  sometimes  unpredictable  ways. 

Second,  while  frozen  ground  is  less  permeable,  it  is  not  impermeable.  The  mathematical 
description  of  transport  of  solutes  through  frozen  ground  is  inherently  more  complex  than 
unfrozen  ground  for  the  following  reasons: 

1.  The  solvent,  water,  is  partitioned  into  two  phases,  liquid  water  and  ice,  that  are  inti¬ 
mately  commingled  in  the  pores  of  the  ground. 

2.  As  the  ice  forms,  the  solutes  (including  the  contaminants)  are  largely  excluded  from 
the  ice,  concentrating  the  remaining  liquid-water  solutions  in  a  thin  film  at  the  colloid 
surfaces.  The  chemical  potentials  of  the  solutes  in  these  solutions,  used  to  estimate  the 
effects  of  solute-solute  and  solute-surface  interactions  are  not  understood  well  enough  to  be 
modeled  accurately. 

3.  In  unfrozen  coarse-grained  soils,  the  movement  of  solutes  is  controlled  by  the  Darcian 
flow  of  water  in  response  to  gravitational  and  pressure  gradients.  In  frozen  ground,  solutes 
may  also  move  appreciably  in  response  to  thermal  and  osmotic  gradients — transport 
mechanisms  that  are  less  well  understood  and  more  difficult  to  parameterize  than  flow  in 
response  to  gravitational  and  pressure  gradients. 
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Darcy’s  law 

Water  flow  in  porous  media  may  be  described  by  Darcy’s  law,  eq  1.  Cold  temperatures 
reduce  water  flows  through  porous  media.  As  has  been  discussed,  cold  temperatures  in¬ 
crease  the  viscosity  of  water.  Freezing  reduces  dramatically  the  permeability  of  soils  and 
ground.  Geochemical  solutions  do  not  freeze  uniformly  at  0°C.  The  equilibrium  freezing 
temperature  of  pore  water  is  a  function  of  pore  geometry,  freezing-point  depressing  effects 
of  solutes,  and  the  charge  behavior  of  the  porous  matrix  (Everett  1961).  The  quantitative 
description  of  these  effects,  which  is  discussed  below,  is  an  active  research  area. 

Soil-water  retention  curve 

The  permeability  of  frozen  porous  media  is  affected  by  its  liquid  water  content.  Accord¬ 
ingly,  a  discussion  of  the  physics  of  unsaturated  porous  media  is  appropriate. 

Matric  potential 

Two  factors  have  been  identified  in  determining  the  relationship  between  matric  poten¬ 
tial  and  water  content  in  unsaturated  soils.  The  first  is  the  pressure  difference  between  the 
soil  atmosphere  and  the  liquid  water  bound  to  soil  particles.  The  second  factor  is  the  work 
required  to  expand  the  liquid  water/air  interface  as  soil-water  content  decreases.  When  a 
stable  interface  is  established,  these  two  factors  are  balanced  and  the  water  content  is  sta¬ 
ble.  For  a  given  matric  potential,  any  environmental  change  that  changes  one  of  these  bal¬ 
anced  factors  relative  to  the  other  can  be  expected  to  cause  the  water  content  to  change. 

Effect  of  temperature 

An  obvious,  and  pertinent,  example  is  temperature,  which  decreases  the  interfacial  tension 
of  water  against  air.  The  effect  of  temperature  on  the  soil- water  retention  curve  may  be  esti¬ 
mated  via 

SHS),®. 

where  0  is  the  volumetric  water  content  of  the  porous  medium  (dimensionless) ,  and  is 
the  matric  potential  of  the  water  in  the  porous  medium  (Pa). 

For  a  given  pore  radius,  the  relationship  between  capillary  pressure  and  interfacial  ten¬ 
sion  is  direct: 

xj/  _  (2y  waCQSft)  (54) 

r 

where  r  is  pore  radius  (m)  and  (|>  is  the  contact  angle  of  the  air/ water  interface  with  the  solid 
(rad).  The  partial  derivative  in  the  denominator  of  eq  54  can  be  calculated  with  the  formula 
recommended  by  eq  9: 

^  =  Yo7r2(l-|J  *[(uv+  v)(T- Ti)- uTi].  (55) 

Accordingly,  one  would  expect  that  the  following  relationship  would  hold  exactly: 

"3'F(r,e) 

ar 

'fr  watt) 

ar 

In  fact  it  does  not,  requiring  the  definition  of  an  empirical  variable  [G(0)]  introduced  by 
Nimmo  and  Miller  (1986)  to  account  for  the  discrepancy: 


=  M  .  (56) 

Y  wa  (0 
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(where  ti  >  t\). 


(57) 


G(0)  = 


y(fe.e) 

TM) 

T  wafe)  i 
Y  wa(h) 


Soil  freezing  characteristic  curve 

Liquid  water  (albeit  generally  very  small  amounts)  exists  in  frozen  soil  at  temperatures 
significantly  below  0°C.  For  water-saturated  frozen  soils,  an  approximation  of  the  liquid- 
water  content  can  be  calculated  from  the  soil-water  retention  curve,  if  known.  The  physics 
of  the  two  situations  are  similar:  in  unsaturated,  unfrozen  soil,  the  liquid-water  content  is 
determined  by  the  balancing  of  the  pressure  difference  across  the  interface  between  the  soil 
water  and  the  soil  atmosphere  with  the  work  necessary  to  deform  this  interface  against  the 
interfacial  tension  of  water  against  air.  The  validity  of  this  hypothetical  mechanism  for  con¬ 
trolling  the  soil-water  retention  curve  is  supported  by  the  observation  that,  for  a  given 
matric  potential,  the  quantity  of  water  held  by  the  soil  decreases  as  the  temperature 
increases;  that  is,  as  interfacial  tension  decreases  with  increasing  temperatures,  less  work  is 
required  to  deform  the  air /water  interface. 

In  frozen,  water-saturated  soils,  there  is  a  water/ice  interface  between  the  liquid  water 
surrounding  the  soil  particles  and  the  ice.  This  situation  is  qualitatively  similar  to  the  water/ 
air  interface  that  exists  in  unsaturated  soils.  A  schematic  displaying  these  parallel  physical 
systems  is  presented  in  Figure  15.  It  is  expected  then  that,  for  a  single  soil  sample  of  unal¬ 
tered  internal  structure,  liquid  water  contents  in  the  frozen,  unsaturated  state  should  be  a 
function  of  the  pressure  gradient  across  the  water/ice  interface.  This  function  should  be 
similar  to  the  empirical  equation  that  describes  the  soil-water  retention  curve. 


Figure  15.  The  similarity  of  capillary  water  of  unsaturated,  unfrozen  soils  and 
unfrozen  water  in  saturated,  frozen  soils. 


Generalized  Clapeyron  equation 

Here  we  follow  closely  the  lucid,  complete  presentation  of  the  chemical-thermodynam¬ 
ic  theory  of  frozen  porous  media  presented  by  Brun  et  al.  (1977).  Consider  three  phases: 
solid  (s),  liquid  (1),  and  gas  (g),  which  are  in  thermal,  but  not  necessarily  in  hydrostatic, 
equilibrium  (i.e.,  Ts  =  T[  =  Tg  =  T,  but  ps*p\±  pg).  Their  development  begins  with  two  sets  of 
classic  equalities:  the  Gibbs-Duhem  equations,  and  the  Kelvin  equation. 

Gibbs-Duhem  equations.  For  each  phase,  the  Gibbs-Duhem  relation  holds: 

o=sf,dr-^dpg+ixfc4if 

B 

0  =  S*ndT-VIUpl+£4cH  (58) 

B 

0  =  S&  dT-  V*  dps  +  X*|dn| 

B 
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where  S^,  S^,  and  5^  are  the  molar  entropies  (J  K-1  mol’1)  of  the  gas,  liquid,  and  solid 
phases,  respectively;  V^.and  are  their  molar  volumes  (m3  mol-1);  x|,  Xg,  and 
Xgare  the  mole  fractions  of  component  B  in  the  three  phases;  and  ji|,  and  are  the 
chemical  potentials  Q  mol-1)  of  component  B  in  the  three  phases. 

Kelvin  equations.  The  Kelvin  equations  give  the  pressure  gradients  across  the  gas/ solid 
and  liquid/solid  interfaces  in  terms  of  the  differential  geometry  of  these  interfaces: 


pg  -  ps  -  ySg 


dVg 


(59) 


pi  -  pS  —  yls 


^sfm 

dVln 


(60) 


where  A|fmand  (m2  mol-1)  are  the  molar  areas  of  the  liquid/solid  and  solid/gas  inter¬ 
faces,  respectively,  and  and  'fs  (N  m-1)  are  the  corresponding  interfacial  tensions.  Capil¬ 
lary  pressures  of  liquid  phases  can  be  estimated  by  applying  these  equalities. 

If  the  solid,  liquid,  and  vapor  phases  are  pure,  then  eq  67  can  be  rewritten: 


0  =  S;8H2OdT”CH2Odpg+dgHg2O 

(61) 

0  =  SS.H2odT-CH2odp'+d^2O 

(62) 

®  =  ^m,H2od7-  Kn?H2OdPS  +  dP-HjO 

(63) 

where  the  superscript  *  indicates  a  pure  phase.  By  taking  the  partial  derivatives  of  eqs  68 
and  69,  dpg  and  dp1  can  be  related  to  dpP  by 


dps  =  dps  +  d  yss 


,dA2 


Sg 


and 


dp1  =  dps  +  d 


(4)- 


(64) 


(65) 


By  subtracting  eq  67  from  eq  69,  and  69  from  68,  then  subtracting  the  latter  difference  from 
the  former,  the  following  relation  is  obtained  after  rearrangement: 


c*s 

°m,H2Q 


c*8 

^m.H2Q 


c*l 

°m,H2Q 


*Vh2Q 


^^mSH20 


V*g 

Vh2o 

vm,H2Q 


V8 

Vh2o 


\/*s 

Vh2o 


^m,H20  “  ^mSH20 ) 


+s,ir)- 


d  T= 


Knl,H2Q 


''Cl.fyO  -  K?,H20 


'  ,.d Asy 
dV*  , 


(66) 


By  noting  that  V^q  »  V£|  H.,0  and  by  assuming  that 

,dAsse 


ysg- 


-  =  0, 


dV* 

the  following  relation  is  obtained: 
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(■^,H20  -  Sm,H2o)dT-  Kl,H2Od^Y1S  j- 
By  defining  liquid-solid  capillary  pressure  (Pa),  p*s,  as 


„ls  _  1  s 
Pc  ~  P  P 


it  can  be  seen  from  eqs  60  and  68  that 


,1s 


dPc  =d|  y 


Therefore,  the  capillary  pressure  of  liquid  water  in  a  frozen  porous  medium  can  be 
calculated  by 

Pcs  T  C*1  _  C*S 

|dpl?  s  p*s  =  J  m-H»°dr  (70) 

0  To  Kn,H20 


P^-]~MHjP-dT  (?1) 

T0  VH20 

where  Tq  (K)  is  the  freezing  point  of  bulk  liquid  phase  and  A1sS^  h2o  is  the  molar  entropy  of 
melting  for  pure  water  (J  K_1  moH). 

Equation  77  is  usually  simplified  by  assuming  that  both  VJ*JH20  and  A]sS*  >H2o  are  con¬ 
stants,  becoming 

p.  =  ^Wt  (72) 

VH20 

Equation  71  ip  often  referred  to  as  the  generalized  Clapeyron  equation.  The  value  of 
A'sS;.H2o/Vm.  Hz0can  be  calculated  from  physical  measurements  reported  in  the  pub¬ 
lished  literature.  Since 

a!  c* _ ^s-^m,H20  (7q\ 


the  molar  enthalpy  of  melting  for  water  measured  byHaida  etal.  (1974),  6006.8JmoH,  can 
be  used.  The  value  of  V£lat  0°C  and  atmospheric  pressure  is  18.0183  cm3  mol-1  (Haar  et  al. 
1984).  The  constant  /  V^Hz0  therefore  is  1.2205  MPa  K-1. 

Following  the  suggestion  of  Koopmans  and  Miller  (1966),  the  Clapeyron  equation  has 
long  been  applied  to  soil  systems  to  estimate  the  matric  potential  of  the  liquid-water  frac¬ 
tion  of  frozen  soils. 

This  argument  begins  with  the  standard  definition  of  capillary  pressure,  pc  (Pa): 

pc  =  p(liquid  soil  water)  -  p(ambient  atmospheric  pressure)  (74) 

which  is  related  to  matric  potential  by 

¥  =  (75) 

PS 

(Scheidegger  1 974) .  A  similar  pressure,  p/,  may  be  defined  for  frozen  ground  to  express  the 
pressure  gradient  across  the  ice/water  interface  in  frozen  ground 
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Pi  =  p(liquid  soil  water)  -  p(ice).  (76) 

A  term  similar  to  matric  potential  can  be  defined  for  the  pressure  gradient  across  the  ice/ 
water  interface  in  frozen  ground: 


(77) 
P  8 

If  it  is  assumed  that  the  volumetric  content  of  soils  for  a  given  matric  or  supercooling  poten¬ 
tial  is  dependent  solely  on  the  ratios  of  pressure  gradients  to  interfacial  tensions, 


then  moisture-release  and  soil-freezing  curves  should  be  related  by 


/ 

eOF)  -  e 

V 


7  wa  Pi  ' 

7 iw  pU> 


(79) 


The  value  of  p2-  as  a  function  of  temperature  can  be  estimated  by  eq  72  and  73  if  it  is  assumed 
that  pj  =  dp1+s  and  t  -  0  =  dT.  The  necessary  data  are 

ywa  (p  =  0.1  MPa,  t  =  0.0PC)  =  0.07564  N  nr1  (Haar  et  al.  1984) 

yiw  (p  =  0.1  MPa,  t  =  0°C)  =  0.033  N  m_1  (Ketcham  and  Hobbs  1996) 


343.2  JmoF1 


-  18.018  cm3  moF1 


V£s=  19.650  cm3  moF1. 

This  implies  that  the  pressure  gradient  across  the  ice/water  interface  changes  with  temper¬ 
ature  according  to 


^!^-  =  7699PaK-i  s7699Pa°C-i. 
dl 

Making  the  proper  substitutions  into  eq  75, 


0(XF)  -  e 


7  wa 


A sMm. 


7 iw  T' j\xs As^rn  Pw^ 


the  following  equation  is  derived: 


(80) 

(81) 


epp)  -  e(ft) 


(82) 


where  f  is  a  constant  approximately  equal  to  -1.80  m  H2O  °C  !. 

Changes  in  hydraulic  conductivity  with  soil-water  content 

The  relative-hydraulic-conductivity  model  proposed  by  van  Genuchten  (1980)  has  been 
adopted  by  many  investigators.  With  some  changes  in  notation,  van  Genuchten’ s  model 
represents  reduced  water  content,  0  (dimensionless),  as  a  function  of  soil-water  matric 
potential  by: 


X-i 
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where  a  (m_1)  and  X  (dimensionless)  are  empirical  parameters.  The  reduced  water  content 
is  defined  by 


__  (Q  (84) 

(e,-er) 


where  0  is  the  volumetric  water  content  (dimensionless),  and  0S  and  0r  are  the  saturated 
and  residual  volumetric  water  contents  (dimensionless).  The  resulting  empirical  equations 
for  volumetric  water  content  are: 


for  unsaturated,  unfrozen  soils 


0  =  0r+(05-0r) 


1 


[OFar+l 


x-\ 

X 


(85) 


and,  by  extension,  for  frozen,  saturated  soils: 

X-l 


0  =  0,+(0,~0r) 


1 


c/tar+i. 


(86) 


Similarly,  the  following  equations  may  be  used  to  calculate  changes  in  relative  permeabil¬ 
ity: 

for  unsaturated,  unfrozen  soils: 


Krd ')  = 


2 


(00|/)x  1 

x-\ 

[l  +  (a\|/)X] 

[l  +  (CX'F)A']  2  A. 


(87) 


and  for  frozen,  saturated  soils: 


(88) 


RESEARCH  PRIORITIES 

The  literature  survey  presented  here  has  revealed  important  gaps  in  scientific  informa¬ 
tion  that  will  impede  the  development  of  contaminant-transport  models  for  cold  regions: 

1.  The  Pitzer  model  provides  a  powerful  method  by  which  to  model  the  thermophysical 
behavior  of  electrolyte  solutions  below  the  freezing  point  of  water.  Unfortunately,  there  are 
few  measurements  with  which  to  apply  this  tool.  Physical-chemical  properties  of  aqueous 
electrolyte  solutions  at  subzero  temperatures:  isobaric  heat  capacities,  molar  volumes  and 
viscosities,  are  needed. 

2.  A  consistent  thermodynamic  treatment  of  the  hydrostatics  of  liquid  water  in  frozen 
porous  media,  a  critical  need,  has  yet  to  be  developed. 
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3.  Aside  from  permeability,  the  pressure-saturation  relation  for  fluids  in  porous  media  is 
the  most  critical  physical  function  for  modeling  the  transport  of  fluids  in  variably  saturated 
media.  It  has  yet  to  be  resolved  entirely  how  temperature  affects  this  relationship.  A  satis¬ 
factory  understanding  of  the  effect  of  temperature  on  soil  retention  curves  for  water  and 
other  liquids  must  be  developed. 

4.  A  co-ordinated  research  effort  to  determine  the  temperature  effects  on  the  physics  of 
insoluble  and  sparingly  soluble  nonaqueous  phase  liquids  in  porous  media  should  be  pur¬ 
sued. 

5.  In  frozen  ground,  the  movement  of  solutes  and  water  in  response  to  osmotic  and 
thermal  gradients  may  be  significant.  Relations  by  which  the  transport  of  water,  solutes, 
and  nonaqueous  phase  liquids  in  response  to  thermal  and  osmotic  gradients  may  be  esti¬ 
mated  reliably  should  be  derived  or  developed. 

CONCLUSION 

A  surprising  number  of  cold-regions  areas  in  the  northern  hemisphere  are  contaminat¬ 
ed.  There  is  a  need  to  incorporate  the  fundamental  physical,  chemical,  and  hydrological 
processes  that  govern  contaminant  hydrology  into  the  numerical  models  that  are  used  to 
assist  in  the  remediation  of  these  sites.  It  is  hoped  this  document  assists  in  this  process. 

The  physical  and  chemical  processes  that  determine  the  fate  of  contaminants  in  cold 
regions  are  just  now  being  understood  comprehensively.  As  with  any  new  area  of  under¬ 
standing,  this  is  a  challenging  process,  but  it  is  also  an  exciting  area  at  the  limits  of  our 
understanding  of  the  dynamic  behavior  of  complex,  heterogeneous  materials. 
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